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; ' Une nbuv'ailA methbde .de calcul dc 1 'interaction onde de 

choc-c6uche limite e8t,prSsentee;'>Le concept de couche limite est uti- 
lise pour decrire 1 'Scoulement A.la parpi, 1 'ficoulement externe est 
approxime, par une relation pression-dgviation (diedre tangent en super- 
sonique linearisS) . , . 


.',v'r Les Equations de couche, .limite sont rSsolues en differences 
finies et la question de. l 'existence et; de 1 'unicite de la solution est 
envisagSe pour le problems direct (pression imposee) ou' inverse (Spais- 
seur de deplacement, coefficient de frottement imposes) . De cette ana* 
lyse ainsi que dcs considerations relatives au couplage visqueux-non 
visqueux, la quantite est choisie comme variable de raccordement . 

I c c 

L'algorithme de couplage presente traite implicitement la 

\ 

condition limite aval necessaire pour farmer correctement le probleme 
de la couche limite en \interaction. L'algorithme utilise une technique 
de linearisation de Newton qui permet une convergence rapide. Un example 
d'application est donne. 
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Ab6tAa.ct ' ■ ‘ . 

, A new m&thod to compute tJiz ^hock mv&-boundoA^ ZayeA tnte.- 
AjCLction pfiuejfvtajd^ The.^boundaAy Zayin^ c.onc^pt tt, lued to de^c/Ubz 
thz vtdcoui ^Zouf neoA th& uxlZZ. and thz outzn. invt&ctd ^Zou) t6 caZcuZa- 
ted by a^.pA<u4>uAz-tuAyu.ng AzZat^n {tangejnt Mzdgz oa ZtnzanU-zzd 6upeA- 
■6on^) {^ofL.&mpZMUJty. ., .,^. . , 


Ihz bouhdoAy Zdy^ zquatiom OAZ 6oZvzd tn ^tyUtz dt^f-QAznzz 
ionm and thz qaz6tton o^ thz zxtstzncz and unicUty o^ thz ioZution -u 
zxanUnzd loA thz cUazcX. {pAZ66uAz ^tmpo6zd) oa tnvo/uz (cUdpZacmznt 
thicknzd^, ^AtctZon coz^^tctznt. . . Zmpo6zd) pAobZzm. fAom thz a^oAzmzn- 
j, .■ tloned amZy^AJi tog^hzA^ulth ot6zp(U-tno^(Ud mtchtng con6tdzAotton6 
. . a poAtccuZoA zhotzz. t& mdz ipA piz. maXzkihd voAtabZz (p^u^6 ) . 


An aZgo/Cithii ts' pAZS'znteii to convzAgz thz matzhzd 6oZutton 
which impZicttzZy tncZudzs thz dowhstAzam condition nzczMi>aAy to cZo&z 
thz intZAacting boandoAy ZayzA pAobZzm. Thz aZgoAithm makz6 luz o^ a 
Nzuiton ZtnzjCLAization tzchfUquz Which t6 &hOwn to pAovidz a ^ai>t convzA- 
gZncz. A typicaZ zxmpZz it -pAztzntzd, 
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coefficient de frotteraent 


fonction d 'amortissement 

fonction de melange amont/aval 

coefficient d* expansion du maillage 

pression externe 

nombre de Prandtl = 

nombre de Prandtl turbulent 

densite externe 

nombre de Reynolds = p U- Xf/y 

00 CO ^ eo 

fonction de sensibilite (couplage) 
vitesse longitudinale externe ou discretisee 
vitesse verticale externe ou discretisee 
coordonnee discretisee 

coordonnee transformee y/e ou discretisee 


h enthalpie statique 

I longueur de melange 

p pression interne 

u vitesse longitudinale interne 

V vitesse verticale interne 

y coordonnees 


Y = io^ction d ' intemittence 

6 epaisseur de couche limite 

epaisseur de deplacement 
- 1/2 

Rg ou viscosite turbulente 

a + 

(1 + £ 


6" 
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Pr 
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rT 


3u 


coefficient de frottement y 

3y 


0 angle, des lignes de courant 

y viscosite dynamique 

0 ) coefficient de relaxation 


ABODE 

ay6c 


coefficients matriciels 



numerotation verticale du maillage 

numerotation longitudinale du maillage 

totale 

charniere 

variable externe 

variable interne ou incompressible 

indice d' iteration 

valeur externe 

valeur a la paroi 

valeur a I'infini 
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I. INTRODUCTION 

Diverses methodes ont ete proposees pour le calcul de 1 ' inte- 
raction onde de choc/couche limite. Les plus anciennes basees sur ba re- 
solution des equations integralcs de couche limite sont encore courammen^ 
utilisees, malgre leur caractere approximatif , du fait de leur rapidite 
{1, 2, 3}. Plus recemnient sont apparues des methodes toujours basees sur 
la notion de couche limite mais resolues cette fois en differences finies 
{^» 5}» puis une solution des equations de Navier Stokes complete est. 
devenue possible grace a I’utilisation de moyens de calculs puissants et 
d ’ algorithmes efficaces {6, 7}. 

La methode de calcul presentee dans cet article a pour objec- 
tif principal d'etudier 1' influence des modeles de turbulence sur la qua- 
lite des previsions. Son application a des ecoulements tridimensionnels 
est aussi envisagee. Ces deux objectifs paraissent pour 1' instant diffi- 
cileraent compatibles avec une resolution directe des equations de Navier- 
Stokes qui serable a la limite des moyens actuels de calcul, au moins pour 
des maillages assez serres. D'autre part, les methodes integrales, dont 
I'efficacite est bien demontree, necessitent un apport semi-empirique cn 
plus de la modelisation de la turbulence et par consequent il est diffi- 
cile de tirer des conclusions sur la modelisation de la turbulence elle- 
meme . La methode choisie est de type couche limite mais resolue par dif- 
ferences finies, le caractere approche resident alors seulement dans 
1 'approximation de couche limite elle-meme. 

L' etude de la modelisation de la turbulence dans 1' interac- 
tion onde de choc/couche limite a fait I'objet d'un travail experimental 
{8} mais le calcul presente est pour I'instant limite a un raodele de 
turbulence algebrique classique, un effort important ayant ete accorde a 
1 'optimisation du calcul de couChe limite et de I'algorithme de couplage. 

Dans le premier paragraphe on rappelle brievement le principe 
de 1 'approche couche limite qui decoule du comportement des equations de 
Navier Stokes aux grands nombres de Reynolds. Les divers modes de raccor- 
dement entre la couche limite et 1 'ecoulement externe sont envisages (a 
la paroi, en 6 , en 6) . 



Le second paragraphe est consacre au calcul de la couche 
limite et en premier lieu au choix du mode de calcul, "direct" ou "inverse", 
et done de la condition limite imposee a la couche limite. Cette etude, 
purement numerique, montre 1' importance du choix de cette condition sur 
1 'existence et I'unicite de la solution locale des equations de coucHe 
limite, en particulier au decollement. 

L' analyse numerique du programme de calcul de couche limite 
est presentee dans le paragraphe suivant. La methode d' integration de la 
region decollee est egalement decrite ainsi que le modele de turbulence 
algebrique utilise. 

Le quatrieme paragraphe est enticrement consacre a 1 'algori- 
thrae de couplage entre la couche limite et 1 'ecoulement exterieur. Une 
methode originale a ete developpee. II a ete demontre {9} que pour etre 
■'bien pose" le probleme de couche limite en interaction forte necessite 
une conaition limite aval. L'algorithme developpe traite implicitement 
cette condition aval et optimise en meme temps le processus iteratif 
de couplage. 

L' application de la methode decrite a un cas d ' interaction 
onde de choc/couche limite en ecoulement supersonique est presentee en 
comparaison avec les resultats experimentaux de LAW {10}. et de divers 
calculs publics a ce jour (5, 6, 7). La prevision des distributions de 
pression parietale et de Cf presente des divergences sensibles avec des 
resultats obtciius par methode de Navier Stokes, en particulier au niveau 
de la pression "plateau". II ne semble pas encore possible d'attribuer 
ces divergences au modele de turbulence utilise, a 1 'amelioration du 
maillage par rapport aux methodes de Navier-Stokes , ou a 1 'approximation 
de couche limite. 

I I . CONDITIONS DE RACCORDEkKNT ENTRE LA COUCHE LIMITE ET L’ ECOULEMENT EXTE KNE 
, II. 1. Generalites 

La methode de calcul presentee dans cet article est de type 
couche limite, e'est-a-dire que I’on admet implicitement que les effets 
visqueux contenus dans les equations de Navier-Stokes sont confines dans 
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une petite region proche de la paroi a I'exterieur. dc laqudle le fluide 
pent etre correctement decrit pair les equations d 'Euler. Cette structure 
des couches visqueuses a des nombres de Reynolds eleves est observee expe- 
rimentalement depuis tres longtemps mais une formulation mathemati que 
rationnelle du comportement des equations de Navier-Stokes Icrsque 
Rg " n'a ete publiee que beaucoup plus tardivement {11, 12}. 

La methode des developpements asyraptotiques raccordes consiste 
a chercher une solution reguliere au voisinage dt e = 0, e = carac- 

terisant les effets visqueux. On recherche done une solution interne 
valide au voisinage de la paroi et une solution externe valide partout 
ailleurs. Les develcppenients asymptotiques externes : 

U(x,y,e) = Ui(x,y) + eU 2 (x,y) + 0(6^) 

V(x,y,e) = Vi(x,y) + eV 2 (x,y) + OCe^) 

R(x,y,e) =■• Ri(x,y) + eR 2 (x,y) + O(e^) 

P(x,y,e) = Pi(x,y) + eP 2 (x,y) + O(e-) 

sont introduits dans les equations de Navier-Stokes et fournissent au 
premier ordre les Squations d 'Euler, et au deuxieme ordre une forme de 
petites perturbations des equations d'Euler. 

Les developpements limites valides au voisinage de y = 0 
sont introduits de la meme maniere apres une transformation de la variable 
normale a la paroi definie par : 

Y = y/e 

y(x,y,e) = Pi(x,Y) + eu 2 (x,Y) + O(e^) 
v(x,y,e) = evi(x,Y) + e2y2-(^.Y) + OCe^) 

P(x,y,c) = Pi(x,Y) + ep 2 (x,Y) + 0(c;2) 

P(x,y^e) = Pi(x,Y) + ep2(x,Y) + O(c^) 

(la forme particuliere du developpement de v est imposee par 1 'equation 
de continuite) . 


( 2 ) 

(3) 


L' injection de ces developpements dans les equations de Navier- 
Stokes fait apparaitre au premier ordre les equations de couche limite de 
Prandtl qui apparaissent done comme une approximation a 1 'ordre 1 des 
equations de Navier-Stokes et au deuxieme ordre une forme de petite per- 
turbation des equations de Prandtl (en 1 'absence de courbure) . 



II. 2. Relations de couplage 


Les conditions de raccordement entre la solution interne et 
la solution externe fournissent les relations ds ccuplage aux divers ordres 
d' approximation {12, 13}. Lorsqu’il s'agit d'une interaction "faible" 
la solution externe d'ordre 1 peut etre calculee en utilisant la condition 
de glissement a la paroi (V(x,0) =0). La solution interne est ensuite 
calculee puis la solution externe d'ordre 2 et ainsi de suite suivant le 
principe de la hierarchie de la solution externe (fig. 1). 


Ce processus tombe en defaut lorsque des chocs apparaissent 
dans la solution externe d'ordre 1 (dP/dx -y “) ou un decollement dans la 
couche limite d'ordre 1 (d^6 /dx^ . Ce type de singularite peut etfe 

evite en considerant simul tenement les etfets d'ordre 1 et 2 a I'exterieur 
avec la couche limite d'ordre 1 (fig. 2). Cette methode de couplage rela- 
tive aux interactions "fortes" fournit la relation de raccordement sur v : 

* ’ •'V* 

V^(x,0) = lim c(vi(x,Y) - Y ~ (x,0)) + (4) 

Y-xo ^ 

ou V designe la solution externe d'ordre 1 et 2 definie par : 

V* = Vi+eV2 


En utilisant 1 'equation de continuite la relation (3) peut etre 
ecrite sous la forme bien connue : 

RV(x,0) = — (ru(x, 0) 6^) R = densite dans I'ecoulement (5) 

externe 

H* 

si I'on cherche a ecrire cette relation a y = 6 par un developpement 
limite au voisinage de y = 0 : 

RV(x,6*) = RV(x,0) + 6* (x,0) + ... 

dy 

et en appliquant I'equation de continuite la relation de couplage s'ecrit : 

’i' 

RV(x,6*) = RU(x,0) ^ 

soit en approximant RU(x,6 ) par RU(x,0) 

RV(x,6*) = RU(x,6*) * ~ 


ou 





( 6 ) 


On peut de la meme maniere obtenir une relation dt. couplage en y = 6 
(14). 
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L'avantage evident du couplage en 6 est que la surface de 
deplaceinent apparait coinme une ligne de courant pour I'ecoulement exte- 
rieur. Par contre daus le couplage en paroi , ,1a couche limite joue le 
role d'une distribution de sources placees a la surface {12} et il y a 
par consequent echange'de fluide entre la couche limite et I'ecoulement 
exterieur (fig. 3). Cette permeabilite de la couche limite pose un 
probleme de raccordement lorsque le fluide externe ne verifie pas la 
condition isentropique valable a 1 'exterieur de la couche limite. Cette 
incoherence ne semble pas cependant remettre en cause la procedure pui^s- 
que des conditions de raccordement supplementaires peuvent etre fournies 
par 1 'etude du raccordement de la couche limite du second ordre (26). 

Une difficult^ associee a la notion de surface de deplacement 
et done au couplage en 6 est le comporteraent supersonique de la couche 
limite au dela d'un nombre de Mach de 1 'ordre de 1.4 pour une couche 
limite turbulente {15}. La couche limite raccordee a I'ecoulement externe 
en 6 ne permet pas alors la naissance du phenomene d ' auto-induction et 
les solutions de branchement correspondant a 1 ' interaction forte n'appa- 
raissent pas. La technique couramment utilisee consiste a appliquer a la 
couche limite une discontinuite (jump supercritique-subcritique) pour deve- 
lopper des solutions d ' interaction forte. Le caractere empirique de cette 
methode peut etre evite en prolongeant I'ecoulement externe a la paroi. 

LE BALLEUR {16} a montre que le couplage en paroi, qui decoule directement 
de la methode des developpements asymptotiques raccordes, garantit, au 
moins dans une plage etendue de nombres de Mach, le comporteraent subcri- 
tique de la couche limite en interaction. Malgre les difficultes liees 
a la permeabilite de la codche limite, ce type de couplage a ete prefere 
pour eviter le recours a la technique du ''jump". 

III. CALCUL DE LA COUCHE LIMITE 
III.l. Introduction 

Parmi les problemes lies au calcul de la couche limite soumise 
a de forts gradients, on peut citer celui de la modelisation de la turbu- 
lence d'une part et d'autre part le calcul du decollement qui implique 


une modification du mode d ' integration, la couche limite n’ayant plus 
dans ce cas le caractere purement parabolique d'une couche limite attachee 

Par opposition aux methodes de calcul dites "directes" dans 
lesquelles le gradient de pression est impose a la couche limite {17} sont 
apparues des methodes "inverses" {18, 19} le gradient de pression etant 
alors le resultat de 1 ' integration des equations avec une condition limite 
externe telle que I'epaisseur de deplacement 6 ou le coefficient de frot-- 
tement parietal Cf. Le choix de la condition externe a prescrire est 
important et la question de I’existence et de I'unicite de la solution 
du probleme de couche limite doit etre envisagee pour les differents 
modes d ' integration . 

III. 2. Choix de la condition limite 


Lors des iterations de couplage entre la couche limite, on 
est amend a imposer des conditions "non physiques" a la couche limite et 
le probleme se pose de savoir si a une abscisse donnee probleme 

de la couche limite a une et une seule solution pour toute valeur de la 
condition imposee. Afin de repondre a cette question la procedure sui- 
vante a ete mise en oeuvre : la couche limite est calculee, jusqu'a 1' ab- 
scisse Xj et au point on fait varier la condition imposde (10). 

L'etat de la couche limite au point J conditionnc en pnrtie le resultat 
et le calcul a ete reproduit pour differenter situations (couche limite 
acceleree, ralentie ou a gradient de pression nul jusqu'a Xj) . Pour 
connaitre le comportement de la solution, decollement inclu^ le calcul 
est effectue soit en methode inverse, (6 prescrit) lorsque I'ecoulement 
est decolle soit en methode directe lorsque I’ecoulement est attache. 

Les termes de convectifs (pu ^ et pu -^) sont negliges lors- 
que la couche limite est decollee afin de preserver la stabilite du calcul 
(Approximation suggeree par Flugge Lotz et Reyhner {21}). Le resultat 
presentc fig. A est relatif a une couchc limite turbulente "supercritique" 
(M = 2.96). La procedure de calcul est la meme que celle decritc au para- 
graphe III. 3. 


En abscisse est porte le rapport p(X )/p(Xj) et en fonc- 
tion de cette variable est representee la variation du Cj^ , de 6 et 
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p u )|:6 , et de p V c'est a dire le flux de masse qui entre ou sort de 
^e e ’ e e 

la couche limite. 

La premiere constatation est 1 'existence d'un rapport de 
pression maximal entre les points J et J+1, le maximum correspondant 
grossierement au decollement de la couche limite . Une fois le decolle- 
ment naissant depasse, la pression diminue legerement lorsque la couche 
limite est plus decollee. Ce resultat est observe quel que soit I’etat 
initial de la couche limite^ qualitativement du raoins^ et montre que : 

1°/ le probleme direct (grad p impose) n'a pas toujours de solution 
2°/ au voisinage du decollement il existe une solution double 

Ainsi le probleme de la singularite au decollement peut etre 
facilement relie a un probleme d 'existence et d'unicite de la solution 
et 1 'utilisation de methodcs directes pour calculer des couches limites 
ddcollees parait assez dangereuse. 

Si I'on considere maintenant le Cj- coinme variable imposee 
a la couche limite on observe que le Cf decroit de. maniere monotone jus- 
qu'au decollement et qu'en ce point la derivee dCf (J+1) /dp(J+l) est infi- 
nie. Dans le doraaine decolle le Cf presente une decroissance monotone 
jusqu'a une valeur limite inferieure et ensuite augmente legerement. 
Notons que rechercher une solution par methode inverse consiste a recher- 
cher 1 ’ intersection de la courbe avec une droite horizontale. On pent en 
conclure que le decollement est facilement f ranchissable a Cf prescrit 
mais que la encore un probleme d 'existence et d'unicite se pose pour de 
forts decollements . 


La variable P^v^ est indiquee mais vQ 1’ allure de la courbe 
PpV^(J+l) = f(p(J+l)) qui presente presque un point de rebroussement , 

1 'emploi de conmie variable est a proscrire absolument. 


* 

ment 6 
6* (J+1) 


Restent les deux variables basees sur I'epaisseur de deplace- 
et le produit . (Fig. 4). La pente negative de la courbe 

= f(p(J+l)) est symptomatique d'un ecoulement supercritique . 


'Ell effet, une couche limite sujiercritique est caracterisee 
par son inaptitude a generer son propre gradient dc pression {2.2,}; Si 



une perturbation externe entraine une augmentation de 1 'epaisseur de 
deplacement soit d6 /dx > 0, le gradient de pression correspondant est 
negatif (grad p<0) . Ceci conduit a un amortissement de la pression ini- 
tiale et la solution couplee ne se branche pas sur des solutions de forte 
interaction {22}. 

NoLons neanmoins qu’au dela d'une certaine compression la 
derivee d6 (J+1) /dp(J+l) redevient positive et les techniques de "jump" 
consistent a passer artif iciellement de la branc’ie supercritique a la 
branche subcritique (11.2). 

En ce qui concerne le calcul de couche limite proprement dit 
il n’a pas eta possible de calculer la couche limite a 6 prescrit dans 
la branche subcritique et de plus un tel calcul serait sujet a caution 
etant donne 1' existence d’une solution double. Par centre dans la branche 
plus verticale (decollement naissant et decollcrnent franc) le calcul a 
6 prescrit ne pose pas de probleme, la solution existe et est uniaue 
mcme pour de tres forts decollements . 


L'utilisation de p u comme variable est plus avantageuse 

car la branche supercritique n'existe pas dp^u^6 (J+1) / dp(J+l) est 

toujours positif et dans le deco) lement la solution est bien definie, 

t 

comme pour 6 . 


Le couplage en paroi faisant intervenir (p^u^6 ) et non 6 
seul^ ainsi que 1 'absence de branche supercritique dans la solution locale 
du probleme de couche limite, conduisent done naturellement a utiliser 
^e'*e^ comme quantite a imposer a la couche limite. 


III. 3. Integration numerique des equations de couche limite 


Le systeme d 'equations a resoudre est le suivant 


a/ 

b/ 

c/ 


'V 


8pu ^ 9pv ^ „ 
9x 9y 


'V/ 


(v = V + 
9u 


p -V 


9u ''' 9u _ ""‘e 9 .V' “ 9u,. 

9x 9y ^e% 9x 9y 9y 


9h 'X- 9h 


9u 


-up u 


e . I- (- " 


. 3h 

e 


e e Dx 9y 


) + 


VI -,9u.?. 
e(-Arr) 


R 


'9y' 
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d/ 

e/ 

f/ 


g/ 


|P = 0 

dy 


(p u -pu)dy = p u 6 = f(x) donne 

e Q. ^ ee 


E = JL. 

P Y-1 

avec les conditions aux limites suivantes 
u(x,0) = v(x,0) = 0 adherence 


3y 


(x, 0) = 0 paroi adiabatique 


u 


he + -^ = h.j. 


enthalpie totale constante a I'exterieun 


( 7 ) 


Le modele de turbulexice e,e est de type algebrique et sera 
discute ulterieurement (III. 5). 


Les equations sont ecrites sous forme adimensionnelle mais en 
variables non transformdes afin de se preter facilement a I'emploi de 
modeles de turbulence plus complexes que les niodeles algebriques. Les 
variables principales choisies sont u et h et les termes nontenant soit p 
soit p sont expriines en fonction de h et hg (enthalpie a I'extcrieur de 
la couche limite) . Les equations sont ensuite completement lincarisecs 
et ecrites sous forme de differences finies du second ordre en x et y. 

Une formulation implicite est evidemment adoptee pour sa stabilite. 


La condition integrale /(p^u^-pu)dy = f(x) est traitee 
simultanement avec les equations de quantite de mouvement selon x et 
d'energie. Contrairement aux. methodes directes la pression et par conse- 
quent sont traitees comme des inconnues. L' integration simultanee des 
equations de quantite de mouvement et d'energie est imposee par la pre- 
sence de p et u dans I'expression de Pg'^g'^ • 

Une integration successive de 1 ’equation de quantite de mou- 
vement et d'energie imposerait une solution iterative de 1’ equation 


a/ 

b/ 


/(PgU^-pu)dy = f(x) du type 

k-1 k k-1 k. , c 

j (Pg u^-p u )dy .= f (x) 


( 8 ) 

ou k designe 1' indice d' iteration 
sur les termec non lineaires 



Cette procedure qui consiste a integier 1' equation de quantite 
de moUvement avec la relation integrale 8b/ ne permet pas le calcul de 
regions decollees a des nombres de Mach eleves car le signe de 


df (x) 


df(x) dp(x,y) 

depend du nombre de Mach local M(x,y) alors que est toujours 

negatif. 


La relation (8. a) est done resolue avec les termes faisant 
intervenir u et p (done h et h^) sous forme explicite et linearisee 

+ + *f 

completement suivant le principe ((AB) = AB+(A -A)B+A(B -B)) 


La structure generale de la matrice pour un profil donne est 
representee fig. 5. Chaque terne represente une matrice Le vecteur 
representant la solution est ecrit sous la forme 


U(l) 

H(l) 

U(2) 

H(2) 


(9) 


U(NDIM) 
^ H(NDIM) 


La matrice est de type tridiagonale avec en plus une ligne et 
une colonne. La derniere ligne est utilisee pour ecrire la relation inte- 
grale 7c/ alors que la colonne represente I'inconnue u(NDIM) = Ug et 
H(NDIM) = hg. 


Chaque ligne est done de la forme 


[A(I)M 


'u(i-i) 

H(I-l) J 


+|B(I) U 


U(I) 


liKDJ 

la derniere ligne s'ecrivant : 
MDIM 

I jE(I)|x 
1^1 


+lc(i)j^ 


U(I+1) 

[Hd+Dj 




u(i) 

H(I) 


= (D(NDIM)) 


( 11 ) 


) 


La resolution est effectuee par un 
derive de I'algorithme de THOMAS {23}. Dans un 
est raise sous la forme presentee fig. 6, ou (l) 
et les coefficients [aj |y| et (6) sont donnes 
substitution : 


algorithme de substitution 

premier temps, la matrice 

, . . . 10 
est la matrice unite ^ ^ 

par les relations de 
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la(I)l = lG(Dl!t:(lB(I)l-iA(I)|*(a(I-l)l)”^ 

|y(I)| = f|G(I)|-|A(I)|>i.•|Y(I-l)l)!t^{|B(I)|-|A(I)|!|:la(I-l)|)■^ (12) . 

(6(I))= ((D(I))-1a(I) |>K6(I-l)))i:(|B(I) l-(A(I))>|r|a(I-l) j)"’' 

La deuxieme ligne est elle-meme reduite par une relation de substitution. 

|e(I)| = |E(I) |-|e(I-l) |t|a(I-l)l 

NDIM-1 

|e(NDIM)| = |E(NDIM)| - J je(I)lt(6(I)) (13) 

1=1 

NDIM-1 

(6(NDIM)) = (D(NDIM)) - I 1 e (I) 1 1 (6 (I) ) 

1=1 

Dans un second temps le vecteur X est calcule par resolution 
de la matrice en commemjant par la derniere ligne et en remontant. Les 
termes non lineaires sont resolus pour chaque profil avant de calculer 
le suivant en iterant jusqu'a convergence. Une sous-relaxation est 
necessaire au moins pour les profils decolles et la valeur relaxee du 
vecteur X a 1' iteration k est donnee par 


Le coefficient de relaxation est choisi egal a 0.5 ce qui cor- 
respond a une optimisation pour 1 ' integration d'une couche limite compor- 
tant des profils decolles. La valeur de u> peut etre prise egale a 0.8 ou 
0.9 si la couche limite est completement attaches. Le test de convergence 
est effectue sur les U(l) . . . .U(NDIM) et 1' iteration arretee lorsque 


U^(I)-U^~^(I) 

U*^(NDIM) 


< ..0001 


pour tous les I. 


Le nombre moyen d' iterations par profil est de I'ordre de 25 

pour le cas de calcul presente. Ce nombre doit pouvoir etre netteraent 

reduit par 1 ’utilisation de modeles de turbulence ne comportant pas de 

termes du genre ou bases sur ce qui n'est pas sans poser des 

problemes lorsque x -K) (Cf. III. 5). 

w 



III. 4. Integration dans la region decollee 


En 1 'absence de decollement 1 ' integration de la couche limite 

dont les equations sont de type parabolique s'effectue dans le sens de x 

croissant, le resultat etant entierement fonction des valeurs initiales 

et des condions aux limites. Le caractere parabolique en +x est perdu 

9 u ^ li 

dans le decollement car les temes pu — et pu — changent de signe ce 

oX 9x 

qui autorise une influence de I'aval vers I'amont. Les equations snnt 
encore de type parabolique dans les regions decollees mais dans le sens 
des X decroissants. 


Trois solutions sont possibles pour- traiter cette influence 
de I'aval vers I'amont : 

1/ 1 'approximation proposee par Flugge Lotz et Reyhner 
2/ le changement du sens d ' integration dans la region decollee 
3/ les balayages successif s de toute la couche limite en modif iaut les 
molecules de differences finies dans les regions decollees 


a) Approximation de Flugge Lotz et Reyhner {21} 


Cette approximation consiste simplement a negliger les termes 
convectifs en x lorsque la vitesse longitudinale devient negative. Ceci 
permet de preserver le caractere parabolique en +x et d'integrer la couche 
limite en un seul balayage en x. Les termes 

3u . 9h . ^ ^ ^ 

pu — et pu — sont transformes en 

oX oX 


p( 

p( 


u+|uL du 


dx 


u+ 1 u L dh 
2 ^ 3x 


(14) 


qui s'annulent done lorsque u<0. Cette methode^ tres utile pour initiali- 
ser une distribution de vitesse dans la couche limite ou pour calculer de 
faibles decollements^ pose neanmoins des probliimes de convergence pour de 
forts decollements . Elle n'a de plus qu'un caractere approximatif , assez 
raisonnable en ecoulement laminaire ou les vitesses de retour sent faibles, 
beaucoup plus contestable en turbulent ou des vitesses de retour de ^0.2 Ue 
sont rencontrees. 
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b) Changement du sens d' integration dans la region decollee 

La solution adoptee par CEBECI {24} d'integrer I'exterieur du bul- 
be de I'amont vers I'aval et le bulbe de I'aval vers I'amont et d'iterer 
jusqu’a convergence du champ global est certainenient la plus efficace sur 
Ife plan du temps de calcul. Cependant I'algorithme est relativement 
complexe. Ue plus dans un calcul couple couche limite/f luide parfait, des 
iterations de couplage sont imposees et cette propriete pcut etre utilisee 
pour converger a la fois le bulbe de decollement et la solution couplee. 

c) Balayages successifs amont-aval 

L'influence de I’aval vers I'amont est dans ce cas rendue 
possible par la modification du type de dif f erenciation ("backward" si 
u(x,y)>ol "forward" si u(x,y)<0) des derivees partielles — et . Pour 
stabiliser la frontiere u(x,y) = 0 une dif fSrenciation centree est appli- 
qude lorsque [ u(x,y)| < u^,. La remontee des perturbations de I'aval vers 
I'amont n'oxcede pas alors 2 mailles par iteration pour un schema du 
second ordre et un nombre d' iterations important doit etre effectue (I'or- 
dre de grandeur est donne par le nombre de profils decolles) . Get inconve- 
nient par rapport a la deuxieme roethode n'est pas critique lorsque des 
iterations de couplage sont a envisager. Le mode de dif ferenciation adopte 
est le suivant : 


3x 


S3 


F+1 ^ 
2 3x 


b 


+ 


F-1 

2 


3u 

3x 


(15) 


avec 

iu ^ 3U(J.I)-4U(J-1,I)4U(J-2,I) 
3x b 2Ax 


^ ^ -3U(K,I)+4U(J+1,I)-U(J+2,I) 

3x f 2Ax 


la fonction F assurant le melange amont aval etant definie par : 


F(J,I) 


(J,I ) 

u(j,iTT 


(l-exp 


( |U(J,I) I) ) 


La vitesse de reference u est egale a ; 

m ^ 


( 16 ) 


u = .02tu 
m e 



Une transition "douce" entre les deux modes de differenciation 
evite une oscillation de la ligne u(x,y) = 0 au cours des iterations de 
convergence des termes non lineaires. La vitesse de reference est choisie 
suffisamment petite pour eviter de prendre en compte des derivees aval 
pour des profils bien attaches. 


III. 5. ilodele de turbulence 

Le modr-le algebrique de turbulence est un modele a deux couche 
classiques. Les termes e et e apparaissant dans 7b et c sont donnes par : 






(17) 


ou e designe la viscosite turbulente, le nombre de Prandtl est pris egal 
a 0.725 et le rapport au nombre de Prandtl turbulent est de Pr/^rt “ 0.80. 

Le modele externe e est donne par : 

0 

e = .0168 u *6? (18) 

0 e 1 

6^ est 1 'epaisseur de deplacement iiacompressible comptee & partir de I'or- 
donnee y ou la vitesse s’annule dans le cas d'un profil decolle 


6* = /^®(1 - f-) dy (19) 

y(u=0) e 

La function d ' intermittence y est representee par 

Y = (l+5.5(y/6)^) 6 etant I'epaisseur do couche limite a 0.99 u^. 

Contrairement aux calculs effectues avec les equations de Navier-Stokes 
completes ou la vitesse varie a 1 'exterieur . de la couche limite {7} et 
ou des chocs prennent naissance au voisinage de y = 6, la frontiere 6 est 
bien definie. 


E . 

1 


Le modele interne est de 


= 1^1 


type longueur de melange 


avec J. = 0.^ y D 

oQ D est le facteur d' amor tissement de Van Driest donne par 1 ’expression 


D . l-e‘5'^5'! 


= epaisseur de la sous couche "laminaire" 
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26y 



Le raodele de turbulence ainsi defini n'a pas fait I'objet de tentative 
d 'optimisation. II est qualifie de modele "d 'equilibre" suivant la 
terminologie de 3HANG et HANKEY et de MAC CORMACK {6, 7}. 


IV. ALGORITHME DE COUPLAGE 

IV . Calcul de I'ecoulement externe 

Bien que le mode de calcul de couche limite soit de type in- 
verse, la relation de couplage permet d'utiliser indif f eremment une methode 
directe ou inverse pour le calcul de I'ecoulement externe. II est prevu 
de coupler le present calcul de couche limite a un calcul externe en me- 
thode des caracteristiques . On utilise actuellement une relation 
pression deviation donnee par : 

P = P (1+A0+B02+C03) (22) 


ou 0 est 1 'angle de I'ecoulement externe prolonge en paroi et deduit de 
la relation de couplage (17) : 


0 " dx 


RU(x,0) dx 


RU(x,0) 


sort 


dx 


dx 


P u 
e e 


(23) 


Ces coefficients dans le developpement limite (22) sont respectivement 


egaux a 
A = 






B = 


GO 

Y+1 


(24) 


C = 5 m2 


IV. 2. Methode de co uplage 

L'algorithine de couplage a pour function de permettre de cal- 
culer la distribution de pression qui satisfait siinultanement les equations 
de la couche limrte et celles de I'ecoulement externe^ les deux ecoulements 
etant relies par la relation de couDlaee (171. 


De nombreuses procedures existent mais les methodes paraissent 
assez laborieuses. Bien que la couche limite soit de nature parabolique 
et que dans le cas d'un ecoulement externe supersonique les equations 
d 'Euler soient hyperboliques , le probleme couple necessite une condition 
a la limite aval pour etre correctement pose {9, 13}. Cette constatation 
rappelle le caractere initialement elliptique dcs equations de Navier- 
Stokes. Le couplage entre les deux regions externe et interne autorise 
en effet la remontee de 1 * information ce qui implique I'existence d'une 
condition a la limite aval pour fermer completement le probleme. 

L' algor ithme propose par VERLE et VATSA {4} qui prend automa- 
tiquement en compte cette condition aval par 1 ' intermediaire de 1 'equation 
de quantite de mouvement ecrite a la paroi constitue une amelioration sen- 
sible par rapport aux methodes de"tir" precedemment utilisees. La methode 
parait mal adaptee au calcul de couches limites decollees (calcul direct 

de la couche limite) et son application a une couche limite supercritique 

^ • • • 

avec couplage en 6 a jete un doute sur la validite des methodes de couche 
limite en interaction {5}. 


La presente methode comme celle de WERLE et VATSA traite impl’i 
citement la condition S la limite aval et optimise le nombre d' iterations 
de couplage. Elle est applicable aux calculs directs ou inverses, et est 
presentee ici sous la formulation '\nixte" (CL inverse, ext direct) . 


Soit k 1' indice d' iteration de couplage. Une estimation du 
gradient de pression en un point dans la couche limite a 1' iteration k+1 
en function de 1' angle 0 de couplage peut etre fournie par la relation : 


iC 

j k+1 , k ^ 9 grad p^ 

grad p = grad p + — j— (0 -0 ) 

90 ^ 

De la meme maniere on peut relier la pression externe p 
en linearisant la relation pression/deviation (22). 

+ AP (o’^'^^-0^) 


k+1 


(25) 
k+1 


a 1 'angle 0 
(26) 


Les angles 0 ainsi que les gradients de pression dans la couche limite 
etant calcules entre les mailles (derivation en deux points) la pression 
externe p doit elle meme etre definie a la demi-maille soit : 


etc . . . . 



21 , 


Les distributions de pression externes et internes devant etre identiques 

k+1 k+1 

a 1' iteration k+1, on doit done verifier p (J) = p (J) . 

Les relations (25) et (26) peuvent etre coinbinees en eliminant 0 ce qui 
donne sous forme discretisee : 


p^'^^(J)-p^'*'^(J-l)-fp*^(J)-p^'"^(J-l)) . p*^^^^(J)+p^'"^(J-l)-^P*'(J)^^P‘'(J-•l)l 


k+1 


k+l,,v . k+1, 


Ax * 


3 grad p 


30 


2 AP 


ou en notant S le rapport 


S = 




2 AP 


(27) 


la precedente relation s'ecrit : 


_k+l/, nx _ p'"(J)*(l-S)-fp^(J)-p^(J-l)]+2>^S^fp^(J)+p’'(J-l)') 
P - S+1 


(28) 


Cette relation permet d’estimer a la maille J-1 la pression p qui doit 
satisfaire les equations de couche limite et de I’ecoulement exterieur. 
Sachant que de plus on impose que la pressj-on aval soit egale a une cer- 

taine valeur P la relation (28) integree de I'aval vers I'amont fournit 

“ , k+1 

une nouvelle distribution de pression p (J=l . . , JMAX) . 


Dans le present calcul la couche limite est calculee en methode 

k+1 

inverse et done, connaissant p (J),il est facile de recalculer 

k+1 

0 a partir d'une des relations (25) ou (26) . La nouvelle distribution 

de PgU est redefinie par 1 'integration de I'amont vers I'aval de la 

« • « 1 

distribution 0 


.!|!k+l, . rt, . , 

p U 0 (x) = p U 0 (x ) + 

e e ' *^6 0 0 


I 6(x)dx 


(29) 


• • • . k+1 

Le profil initial est perturbe sa pression devant etre egale a p (1) 

mais le produit sst maintenu constant. On retrouve ainsi la solution 

de faible interaction, dans la mesure ou la forme du profil initial est 

compatible avec cette solution. En tout etat de cause un profil de type 

plaque plane fournit une bonne approximation. 


La presente methode se revele extremement rapide pour les 
interactions ne conduisant pas au decollement de la couche limite (15 a 
20 iterations) (fig. 7). Une sous relaxation est necessaire pour des calculs 


incluant des decollements du fait d’une part de la non linearite des 

'i' 

courbes /p au voisinage du decollement (fig. 4) et d'autre part 

parce que le bulbe de decollement est integrc par des balayages successifs 
et que done les distributions de p^u^6 ne doivent pas varier trop d’une 
iteration a 1' autre. Le nombre d' iteration necessaire est alors de 30 a 40. 

IV. 3. Calcul de la derivee S = -■ ^ 

do 

La lineatisation de la relation du comportement local de la 
couche limite utilisee pour 1 'acceleration de la convergence de la methode 
de couplage necessite la connaissance du terme S(x) = — ^ (Cf rela- 

tion 28 ) . 


Notons que pour une couche limite de caractere stricteraent 
parabolique (attachee) la solution en un point donne x = Xj depend, outre 
des conditions initiales, de la distribution cn x de la condition a la 


limite, soit en methode inverse de p u 6^ de x “ X, & x = X, ; 
’ ’^e e 1 J 


p(Xj) “ (p^u^6^(Xj^), L = 1 J) 


(30) 


L'optimisation idSale du processus iteratif de couplage imposerait la 
connaissance des derivees (Jacobien) 
dp, 


dp u 6? 
e e L 


pour L = 1, , . . . , J 


(31) 


ce qui conduirait au calcul de termes. La linearisation de la relation 
grad p(J)/e(J) revient done a ne considerer que les deux diagonales princi- 
pales de la relation (31). Si k est 1' indice de 1' iteration de couplage 
la pression en J a 1' iteration k+1 est estimee par : 


= p’"‘"\j-i) + s^(j)(e^''^(j) -e*"(J)) 


(32) 


On observe pratiquement que la convergence est assuree sans 
relaxation sur la distribution de 6(J) done de (J) lorsque la couche 

limite est attachee. Par contre lorsque le decollement se produit la 
derivee S(J) change de signe dans les regions decollees et une sous relaxa 
tion parait alors necessaire. On doit en plus considerer que 1' influence 
de I’aval sur I’amont remet, au moins en theorie, la validite de la rela- 
tion (31) en cause puisque la pression en J depend alors aussi de la 
distribution de pour L>J. 
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Cependant pour Ic cas du calcul preseute en (V) les derivees 
grad p (Xj) et 9(Xj) ont ete evaluees en negligeant les points en aval, 
soit : 

grad Pj_^/2 

e. 


p(Xj) - P(Xj_l^ 


Ax 




La methode la plus directe pour obtenir la derivee d grad p/d0 
consists a faire deux calculs successifs au point Xj en perturbant la valeur 
de au point J. Cette methode necessite alors deux calculs de c'ouche 

limite par iteration de couplage . 


L 'etude de la relation de comportement (III. 2) a montre que le 
teme de sensibilite S (pente de la normale a la courbe /p) changeait 

de signe lors du decollement. II est alors interessant de relier S au C£ 
pour eviter le double calcul de couche limite : 

s‘"(J)o: C^(J) (33) 

Le coefficient de proportionnalite est calculi en divers points 
du champ par un double calcul de la couche limite comme indique precedem- 
ment et la relation (33) fournit une bonne approximation pour les autres 
points . 


V. RESULTATS ET DISCUSSIONS 


V.l. Resultats 


Le calcul presente est relatif au cas experimental etudie par 
LAW {10} et qui a servi de standard pour divers calculs d' interaction soit 
en methode de couche limite (5, 22} soit par resolution des equations de 
Navier-Stokes (6, 7}. Lemaillage utilise est de 60x60, le pas en x est 
constant, le pas en y variable avec un coefficient d'expansion de K = 1.120 

AY(I+1) = Ki;AY(I) 

La question du maillage au voisinage de la charniere n'a pas 
ete abordee. Le calcul est effectue sur paroi plane, la deviation inter- 
vient simplement au niveau de I’ecoulement externe, c’est-a-dire quo le 
cas de calcul est relatif a 1 'interaction par choc incident et non pas 



par deviation de la paroi. L’ equivalence entre les deux modes de compa- 
raison est bien etablie en laminaire car 1 'interaction s'etend sur un 
domaine relativement grand par rapport a I'epaisseur de couche limite et 
done par rapport a la maille en X. Dans le cas present 1 'equivalence n'est 
pas assuree, comme en temoigne le calcul de SHANG et HANKEY {27} effectue 
pour les deux modes de compression. 


Les donnees du calcul sont = 2.96, le nombre de Reynolds a 

la charniere de Rg = 10^, le nombre de Prandtl Pj- = 0,725. Le domaine de 

calcul s'etend de X /X =0,88 h X /X = 1,124 et Y /X = 0,07. 

o c ’ max c max c 


La viscosite dynamique est approximee par une loi en puissance : 




0.76 


(34) 


Les deux profils initiaux sont donnes par un calcul de couche 
limite direct {25}. Le calcul d ' interaction est initialise par une pre- 
miere iteration effectude en methode directe a gradient de pression nul . 

La pression itnposee dans les iterations de couplage est imposee en 
ct egale a la pression a 1' inf ini aval fournie par la relation (22). Compte 
tenu de la sous relaxation appliquSe dans la relation de couplage la pres- 
sion infini aval n'est atteinte qu'au bout de 7 iterations. La valour impos§e , 
est ensuite depassee puis converge lentement (fig. 8). L'evolution de p et du C 
avec le nombre d'itefations est presentee sur la figure suivante (fig. 9). 

La pression parietale finale est presentee sur la figure 10 en comparaison 
avec les resultats experimentaux et ceux obtenus par divers auteurs {5,6,7}. 

Le Cf final est presente figure 11. 


Le temps de calcul necessaire pour obtenir une convergence de 
0,5% sur la pression finale est de 12 minutes sur IBM 370/168 ce qui cor 
respond a 40 iterations de couplage environ. 


V.2. Discussion 


L'observation des resultats que ce soit la pression parietale 
ou le Cf (fig. 10, 11) montre des ecarts importants entre les divers resul 
tats obtenus. L' approximation de couche limite, dont la validite semblait 
douteuse a la suite des travaux de WERLE et BERTKE, ne semble pas devoir 
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etre remise en cause. II semble que dans le cas cite I’algorithme de cou- 
plage (methode directe/couplage en 6 ) soit plus vraisemblablement a I'ori- 
gine du resultat tres decevant. 

Le modele de turbulence employe dans le present calcul se 
rapproche tres nettement de celui de HUNG et MAC CORMACK, la difference se 
situant au niveau du modele interne ; 

p£^|”| present calcul 

1/2 

HUNG et MAC CORMACK {7} 

'•ay 3x ■' 

Cette difference semble assez mineure et les differences entre 
les calcul de HUNG et MAC CORMACK et le present calcul peuvent avoir plu- 
sieurs origines : 

1) la cause la plus probable est liee a la resolution numerique faite sur 
paroi plane, c'est-a-dire que le cas traite est celui de 1 ' interaction 
par choc incident. Le calcul comparatif d'une interaction par choc in- 
cident et par deviation de la paroi {27} montre que dans le premier cas 
la distribution de pression parietale est decalee vers I’amont, ce qui 
correspond bien au resultat, et que le minimum du C^ avant le recolle- 
ment est lui moins prononce ce qui est egalement observe. 

On peut aussi citer ; 

2) des questions relatives au raaillage (60 pas en y dans le present calcul, 
31 seulement dans le calcul de HUNG et MAC CORMACK pour la couche dissi- 
pative et une part de I'ecoulement non visqueux) . 

3) La definition des grandeurs telle qua u^.la valeur a la frontiere de la 

couche limite, en methode de Navier-Stokes, (SHANG et HANKEY utilisent 

u qui semble se situer en dehors de la couche limite) 
max • ' 

4) L 'approximation dp/dy = 0 de couche limite 

L'obtention des resultats prescntes est assez rScente et il 
n'a pas encore ete effectue de tests comparatifs en particulier 1 'evalua- 
tion des termes negliges dans 1 'equation de quaiitite de raouveraent normale 




a la paroi. Ces controles ainsi qu'un effort sur la modelisation des terraes 
turbulents seront entrepris ulterieurement . La rapidite du calcul autorise 
en effet une etude plus systeraatique des modeles de turbulence dont 1' impor- 
tance est evidente pour le traitement des ecoulements decolles. 


VI. CONCLUSION 

Une methode de calcul de 1 * interaction onde de choc/couche 
limite a dte developpee. Elle utilise le concept de couche limite formulee 
en methode inverse, la condition a la limite imposee est le produit 
Cette condition limite semble constituer le meilleur choix, malgre la rela- 
tive complexite qu'elle entraine au niveau de 1 ’ integration de la couche 
limite, car d'une part le couplage conserve un caractere subcritique a des 
nombres de Mach eleves en ecoulement turbulent, et d'autre part I'existence 
et I'unicite de la relation parait assuree dans les regions de forts gra- 
dients, decollement inclus. 

L'algorithme de couplage base sur le comportement local de la 
couche limite prend implicitement en compte une condition imposee a 1 ’aval- 
pour le probleme coupld. Les durees de calcul sont extremement competitives j 
par rapport a une resolution des equations de Navier-Stokes completes ce qui 
autorise 1 ’emploi d'un maillage serre (60x60) sur un mini-calculateur . 

Des divergences apparaissent par rapport aux methodes exactes (NS) 
mais il n’est pas possible d'en tirer de conclusion definitive, le mode de 
production de 1 'interaction OC/CL n'etant pas le meme (choc incident/devia- 
tion de la paroi) . 

Plusieurs applications sont prevues ou en cours actuellement . 

En premier lieu le raccordement a une methode exacte de 1’ ecoulement externe 
(caracteristiques) . Une extension a des ecoulements tridimensionnels est ega- 
lemcnt prevue. Enfin un gros effort doit etre consenti envers la modelisation 
de la turbulence dont 1’ importance parait primordiale pour une prevision 
"raisonnable" de 1 'ecoulement . 


Ce travail a ete soutenu par la Direction des Recherches, Etudes Techniques 
sous le contrat n° 78/240. 
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FIG.l - Couplage regulier (hierarchique) entre la 
couche limite et l*ecoulement exteme - 



ler ordre 


2§me ordre 


FIG. 2 - Couplage singulier (simultane) entre la 
couche limite et I'ecoulement externe - 
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Relation locale de comportement entre 
P V , 6 , P„u 6 , le C- et la pression 
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Structure de la matrice en methode 
inverse (6 ou p u 6 impose) - 







FIG. 6 - Structure de la 
matrice apres 
reduction 
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FIG. 7 - Convergence de la 
pression finale 
(couche limite 
attachee) 


Nombrc d' iterations 
de couplagc 



FIG. 8 - Convergence de la 
pression finale 
(couche limite 
decollee) 


+ 


10 


20 


30 


Nombre d' iterations 
de couplage 









FIG. 10 - Coniparaison entre les distributions de pression 

parietalc (Modele de turbulence en "equiiibre") “ 




